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The  combination  of  iterative  methods  with  preconditionings  based  on  incomplete  LU  factorizations 
constitutes  an  effective  class  of  methods  for  solving  the  sparse  linear  systems  arising  from  the 
discretization  of  elliptic  partial  differential  equations.  In  this  paper,  we  show  that  there  are  some 
settings  in  which  the  incomplete  LU  preconditioners  are  not  effective,  and  we  demonstrate  that 
their  poor  performance  is  due  to  numerical  instability.  Our  analysis  consists  of  an  analytic  and 
numerical  study  of  a  sample  two-dimensional  non-self-adjoint  elliptic  problem  discretized  by  several 
finite  difference  schemes. 
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1.  Introduction 


The  preconditioned  conjugate  gradient  method  [3,  4]  and  preconditioned  iterative  methods  for 
nonsymmetric  linear  systems  [1,6,  17,  22]  are  effective  methods  for  solving  the  large  sparse  linear 
systems  arising  from  the  discretization  of  elliptic  partial  differential  equations.  Two  good  precondi¬ 
tioners  are  the  incomplete  LU  factorization  (ILU)  [15]  and  the  modified  incomplete  LU  factorization 
(MILU)  [5,  11],  each  of  which  makes  use  of  an  approximate  factorization  of  the  coefficient  matrix 
into  the  product  of  a  sparse  lower  triangular  matrix,  L,  and  a  sparse  upper  triangular  matrix, 
U.  For  the  symmetric  positive-definite  systems  derived  from  the  finite  difference  discretization  of 
self-adjoint  elliptic  problems  on  a  uniform  n  x  n  grid,  it  is  known  that  the  MILU  preconditioning 
produces  a  reduction  of  the  condition  number  from  0(n2)  to  O(n).  The  ILU  preconditioning  does 
not  improve  the  conditioning  in  this  way,  but  it  has  been  observed  empirically  to  generate  a  linear 
system  most  of  whose  eigenvalues  are  clustered  near  one.  The  effectiveness  of  both  techniques  for 
the  nonsymmetric  linear  systems  derived  from  non-self-adjoint  elliptic  problems  has  been  demon¬ 
strated  in  many  numerical  experiments  [2,  7,  8,  20],  although  the  analysis  from  the  symmetric  case 
has  not  been  generalized. 

Let  A  denote  the  coefficient  matrix,  let  Q  =  LU  denote  the  approximate  factorization  of  A, 
and  let  R  =  Q  -  A.  Loosely  speaking,  the  analysis  for  symmetric  discretized  elliptic  equations 
examines  the  effect  of  R  on  vectors  u  whose  values  come  from  a  smooth  function  evaluated  at  the 
mesh  points.  In  particular,  a  heuristic  explanation  of  the  difference  between  the  MILU  and  ILU 
factorizations  is  that  the  individual  entries  of  the  vector  Ru  satisfy  [i?u]j  =  O(lfn)  for  the  MILU 
factorization,  whereas  [.Ru],-  =  0(1)  for  the  ILU  factorization  [11].  The  MILU  factorization  has  a 
higher  order  of  accuracy  as  an  approximation  to  A  than  the  ILU  factorization  (see  also  [18]).  In  this 
sense,  the  analysis  is  reminiscent  of  the  notion  of  “consistency”  of  difference  schemes  for  ordinary 
differential  equations  [10].  The  notion  of  order  of  accuracy  also  extends  to  the  nonsymmetric  case 
(see  e.g.  [14]).  In  this  regime,  the  MILU  factorization  is  also  of  higher  order  of  accuracy,  and  it 
has  been  demonstrated  to  be  more  effective  in  many  numerical  experiments  [7,  8]. 

In  this  paper,  we  show  that  stability  can  also  play  a  role  in  the  performance  of  the  incomplete 
LU  preconditioners  when  they  are  applied  to  discretized  non-self-adjoint  elliptic  equations.  We  show 
using  a  model  problem  that  there  are  nonsymmetric  linear  systems  that  can  cause  difficulty  for 
either  the  ILU  preconditioning  or  the  MILU  preconditioning,  and  that  the  source  of  this  difficulty  is 
instability  of  the  computations  involving  the  triangular  factors.  Our  analysis  is  similar  to  stability 
analysis  for  methods  for  ordinary  differential  equations  [10].  It  shows  that  the  performance  of 


incomplete  factorizations  is  sensitive  to  both  the  values  of  the  coefficients  of  the  elliptic  operator 
and  the  choice  of  difference  scheme  used  to  discretize  the  problem.  In  Section  2  we  present  the  model 
problem  and  give  examples  of  numerical  difficulties  exhibited  by  both  preconditioners.  In  Section  3 
we  construct  constant  coefficient  approximations  to  the  two  factorizations  based  on  an  asymptotic 
analysis  of  the  values  of  their  coefficients,  and  in  Section  4  we  present  the  stability  analysis  for 
these  simplified  factorizations.  In  Section  5,  we  demonstrate  with  numerical  experiments  that  the 
presence  of  instability  correlates  with  a  degradation  of  the  performance  of  the  preconditioners. 


2.  The  Model  Problem  and  Some  Numerical  Examples 

In  this  section,  we  present  a  model  problem  and  briefly  discuss  some  numerical  experiments 
that  demonstrate  some  difficulties  encountered  by  the  ILU  and  MILU  preconditioners.  Consider 
the  constant  coefficient  elliptic  equation 

-Aw  +  2P1uI  4-  2Piuy  =  /  (2.1) 

on  the  unit  square  f2  =  {(or,  y)|0  <  x,  y  <  1},  with  Dirichlet  boundary  conditions  u  =  f  on  dQ. 
Discretizing  (2.1)  on  an  n  x  n  grid  gives  rise  to  a  sparse  linear  system  of  equations 

Au  =  g  (2.2) 

of  order  N  =  n2.  We  consider  two  difference  schemes  for  approximating  the  first  derivatives  in 
(2.1).  In  the  first  scheme,  we  use  second  order  centered  differences 

„  _  ui+l,j  “  wi  — I,/  .  H«./+l  —  Uj  j—i 

*’* - 2 h - '  - 2 h - ' 

where  h  =  l/(n  +1).  In  the  second  scheme  we  use  first  order  upwind  differences,  i.e.  backward 
differences  if  the  coefficient  is  positive,  forward  differences  if  the  coefficient  is  negative.  For  the 
upwind  scheme,  we  restrict  our  attention  to  the  case  Pi,  Pi  >  0,  so  that  the  differences  are  given 
by 

j  —  Uj  —  ij  Ujj  —  Uj  j  —  l 

Uz  — - ; - Uj,  ^  — - - - - - . 

h  h 

For  both  schemes,  we  use  standard  second  order  centered  differences  for  the  Laplacian  [9] 

A  ui+i  j  ~  2 w,;  +  Uj—ij  Wi,;  +  i  -  2utJ  + 

h2  h 2 

Unless  otherwise  specified,  we  order  the  unknowns  along  lines  in  the  x-direction,  i.e.  the  indices 
{ij}  are  ordered  so  that  for  each  j,  contiguous  unknowns  correspond  to  points  whose  i-indices  vary 
from  1  to  n.  The  coefficient  matrix  has  the  form 


a  d 
bad 


bad 
b  a 

c  ad 

bad 

b 


V 

After  scaling  the  matrix  and  right  hand  side  by  h2,  the  matrix  entries  are  given  by 

a  =  4,  6=-(l  +  pi),  c=-(l+p2), 

d  =  —  1  +  pi ,  e  =  - 1  +  p2 

for  the  centered  difference  scheme,  where  pi  =  P^h,  p2  =  P2/z,  and  by 

a  =  4  +  2(pi  4-  p2)  b  =  -(1  +  2pi),  c  =  —(1  +  2p2), 
d  =  — 1,  e  =  — 1 
for  the  upwind  scheme. 

We  note  that  in  practice  the  choice  of  mesh  size  h  may  be  affected  by  accuracy  considerations 
and  the  sizes  of  the  coefficients  Pi  and  P2.  In  particular,  if  the  solution  contains  a  boundary  layer, 
then  with  too  coarse  a  mesh  the  centered  difference  scheme  will  produce  an  oscillatory  solution 
and  the  upwind  scheme  may  not  resolve  the  boundary  layer.  See  e.g.  [19]  for  a  discussion  of  these 
issues.  Our  focus  in  this  paper  is  on  the  algebraic  properties  of  the  matrix  equation  rather  than 
the  accuracy  of  the  discrete  solution.  In  making  this  study,  we  will  have  occasion  to  consider  some 
large  values  of  P,7i  that  may  not  practical  in  applications.  (See  Section  5.) 

We  will  give  precise  definitions  of  the  two  incomplete  factorizations  under  consideration  in 
Section  3.  .4s  evidence  that  their  performance  on  similar  looking  problems  can  vary  greatly,  we 
consider  here  three  instances  of  equation  (2.1): 

Problem  1  :  Pi  =  0,  P2  =  50, 

Problem  2  :  Pi  =  P2  =  50, 

Problem  3  :  P\  =  -50,  P2  =  50. 

In  each  case,  the  right  hand  side  /  of  (2.1)  is  determined  so  that  the  solution  is 

u(x,y)  =  xezy  sin(Trx)sin(iry).  (2.4) 
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bad 
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b  a  J 


(2.3) 
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We  discretize  on  a  31  x  31  grid  using  centered  differences  for  the  first  derivatives.  The  precondition¬ 
ing  is  applied  from  the  right,  i.e.  given  the  incomplete  factorization  Q  =  LUy  the  preconditioned 
problem  to  be  solved  is 

AQ~1u  =  g,  u  =  Q~l  u.  (2.5) 


Using  Orthomin(l)  [6]  as  the  basic  iterative  method,  we  attempt  to  solve  each  of  the  three  problems 
with  both  the  ILU  and  the  MILU  preconditioning.  Table  1  contains  the  number  of  iterations 
of  Orthomin(l)  needed  to  satisfy  the  stopping  criterion  of  ||r,-||/|ro||  <  10~6  (where  |ju||  denotes 
the  Euclidean  norm  (v,  u)1/2).  The  symbol  oo  indicates  that  the  residual  norms  {||r;|||  stopped 
decreasing  at  some  point  of  the  iteration,  i.e.  that  the  iteration  failed  to  converge  to  the  solution. 


ILU 

MILU 

Problem  1. 

21 

18 

Problem  2. 

oo 

7 

Problem  3. 

32 

OO 

Table  1:  Number  of  iterations  of  Orthomin(l)  to  convergence. 


ILU 

MILU 

Problem  1. 
Problem  2. 
Problem  3. 

Leftmost  Rightmost 

0.172,0.196  1.31,1.32 

-49.3.-12.6  18.2,56.1 

0.043.  0.075  1.56,  1.57 

Leftmost  Rightmost 

0.240,  0.247  2.589,  2.593 

0.681,  0.681  1.02,  1.03 

-2.90E10,  -4.36E8  4.36E8,  2.90E10 

Table  2:  Leftmost  and  rightmost  eigenvalues  of  the  symmetric 
parts  of  the  coefficient  matrices. 


Note  that  if  the  symmetric  part,  (AQ~l  +  AQ~l)/2,  of  the  coefficient  matrix  in  (2.5)  is  positive- 
definite,  then  Orthomin(l)  is  guaranteed  to  generate  a  sequence  of  iterates  whose  residual  norms 
are  strictly  decreasing  [6,  7).  The  symmetric  part  is  indefinite  if  and  only  if  it  has  both  nonpositive 
and  nonnegative  eigenvalues.  In  Table  2.  we  list  the  two  leftmost  and  two  rightmost  eigenvalues  of 
the  symmetric  parts  of  the  six  coefficient  matrices  tested.  This  data  confirms  that  the  problems 
in  which  the  failures  occurred  have  indefinite  symmetric  parts.  Taken  together,  the  results  from 
the  two  tables  also  indicate  that,  at  least  when  centered  differences  are  used  for  the  discretization 
of  (2.1),  the  performance  of  the  incomplete  factorization  preconditionings  is  very  sensitive  to  the 
values  of  the  coefficients  of  the  first  derivatives.  (We  remark  that  the  symmetric  part  of  the  original 
matrix  A  is  the  discrete  Laplacian.  which  is  positive  definite  [9j.) 
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3.  The  Incomplete  Factorizations 

In  this  section,  we  define  the  ILU  and  MILU  factorizations  and  construct  simplified  constant 
coefficient  approximations  to  each  of  them  that  will  lend  themselves  to  a  stability  analysis.  We 
consider  incomplete  factorizations  in  which  the  lower  and  upper  triangular  factors,  L  and  U,  have 
the  same  sparsity  structure  as  the  lower  and  upper  parts  of  A,  respectively,  and  U  is  unit  diagonal. 
For  the  five-diagonal  matrix  of  (2.3),  the  factors  have  the  form 


/  <*i 


L  = 


02 

OC2 

03 

'ln+1 

^r»4-2 

/l  8\  hi 

1  <52  T]  2 


U  = 


VN-n 


<!>W-1 

1  J 


&N- 1  I  I  1 

V  7/V  0N  aN  ' 

The  ILU  factorization  is  defined  so  that  the  entries  in  the  product  matrix  LU  are  equal  to 
entries  of  A  whenever  the  latter  are  nonzero  [15] .  Multiplying  L  and  U  and  setting  the  coefficients 
of  the  product  equal  to  the  appropriate  nonzero  entries  of  A  gives  the  following  defining  equations 
for  the  nonzero  entries  of  L  and  U :  . 


0)  ~ 


'  a 

]  =  1 

a  —  0j  8j  _  j , 

2  <  j  <  n 

a  -  IjTlj-n , 

n  +  1  <  j  <  N  and  j  =  1  mod  n 

■  a  ~  0j6i- 1  “  Ijlj-n, 

n  +  1  <  j  <  N  and  j  ^  1  mod  n 

(  b  j  ^  1  mod  n 

(  c  n  +  1  <  j  <  N 

1)  =  i 

( 0  otherwise 

(  0  otherwise 

(3.1) 


*> 


d/otj  j  i1  0  mod  n 

0  otherwise 


1j 


e/otj  1  <  j  <  iV  —  n 


0  otherwise 


The  MILU  factorization  is  defined  so  that  the  nonzero  off-diagonal  entries  of  A  are  equal  to  the 
corresponding  entries  of  LU ,  and  the  sum  of  the  entries  of  each  row  of  R  =  LU  -  A  is  rh 2.  r  >  0. 
We  restrict  our  attention  to  r  =  0.  (This  has  been  observed  to  be  a  good  choice  in  practice  [3.  7. 
12],  although  the  analysis  of  MILU  is  only  known  to  hold  for  r  >  0  [5,  11].)  The  zero  row-sum 
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condition  is  imposed  for  the  j’th  row  by  subtracting  the  error  occurring  in  row  j  of  R  from  the 
j' th  diagonal,  ay.  The  defining  equations  are  as  in  (3.1),  except  the  diagonal  entries  are  given  by 


aj  =  < 


a , 

a  -  pj&j-i  -  Pjrij-i , 
a  -  -  ljm-n  - 

a  —  —  7j >] j-n  ~ 

a  -  an  j-n  -  Ij&j-n, 
a  (3j8j— i  —  ')jrlj—ni 


j  =  1 

2  <  j  <  n 

n+l<j<N  —  n  and  j  =  0  mod  n 
N  -  n  <  j  <  N  (3.2) 

n  +  1  <  j  <  N  and  j  =  1  mod  n 
j  =  N 


[a  -  0j6j-i  -  l,rij-n  ~  0jVj- 1  ~  IjSj-n,  otherwise. 

For  the  ILU  factorization,  all  but  2n  —  1  of  the  diagonals  {<*>}  satisfy  the  last  recurrence  of 
(3.1).  For  the  MILU  factorization,  all  but  4n  —  4  diagonals  satisfy  the  last  recurrence  of  (3.2).*  To 
facilitate  an  analysis,  we  construct  constant  coefficient  ILU  and  MILU  factorizations  based  on  an 
asymptotic  analysis  of  these  two  recurrences.  We  make  use  of  a  relationship  between  the  recurrences 
for  {ay}  and  continued  fractions. 

.4s  motivation,  consider  the  first  block  (i.e.  the  first  n  rows)  of  the  ILU  factorization,  whose 
diagonal  entries  satisfy 


ai  =  a,  ay  =  a  —  bd/aj-i,  2  <  j  <  n. 


(3.3) 


Expanding  for  the  y’th  value  gives 


bd 


a,  =  a 


bd 


>  j-1  divides. 


(3.4) 


’ 1  •  a  —  bd/a  ) 

In  the  language  of  continued  fractions,  is  the  (j  -  l)’st  approximant  of  the  continued  fraction 
([21]  p.  14) 

a  - 


bd 


bd 


a  — 


a  — 


The  convergence  of  continued  fractions  (i.e.  of  the  approximants)  of  this  form  is  well  understood. 
We  state  without  proof  the  result  we  need,  which  is  taken  essentially  verbatim  from  [21]  (Theorem 
3.2,  p.  39). 


'Actually,  with  the  convention  that  the  off-diagonal  entries  are  zero  in  the  appropriate  indices,  all  the  diagonal  entries  satisfy 
the  last  recurrences  of  (3.1)  or  (3.2).  We  specify  the  exceptions  explicitly  to  emphasize  that  there  are  exceptions. 


Theorem  3.1.  The  continued  fraction 


1  + 

converges  for  any  (complex)  number  £  except  when  £  =  —  ~  —  s,  where  s  is  real  and  positive.  For 
real  its  value  is 

1  +  y/1  +  4g 
2 

To  use  Theorem  3.1  to  examine  the  convergence  of  (3.3)  or  (3.4),  consider  ay  =  ay/a.  These 


quantities  satisfy 


bd/a2 

dq  =  1,  a.  =  1 - — ,  2  <  j  <  n. 

<*j- 1 


With  £  =  —bd/a2,  the  approximant  for  ay  analogous  to  (3.4)  is  the  (j  —  l)’st  approximant  of  a 
continued  fraction  of  the  form  (3.5).  The  convergence  result  requires  that  £  >  -1,  which  follows 
for  both  the  centered  and  upwind  difference  schemes  by  direct  computation.  Hence,  {ay}  converges 
to 

1  +  \J  1  -  Abd/a 2 


and  {ay}  converges  to 


wt1)  =s  aa  = 


a  +  Vr  - 


For  the  centered  difference  scheme,  the  limiting  value  is  2  +  \JZ  +  p\,  and  for  the  upwind  scheme 
it  is  2  +  pi  +  P2  +  \/(2  +  pi  +  P2)2  -  (1  +  2pi).  We  remark  that  the  limiting  value  (3.6)  is  equal  to 
the  larger  root  obtained  by  formally  substituting  a  constant  value  in  place  of  a y  and  ay_i  in 
(3.3)  and  solving  the  resulting  quadratic  equation  for  a^1). 

This  argument  establishes  the  convergence  of  the  sequence  defined  by  (3.3)  as  j  —  00,  although 
it  does  not  prove  that  the  first  n  values  are  near  the  limiting  value.  Moreover,  most  of  the  diagonal 
entries  of  L  outside  of  the  first  block  satisfy  the  more  complicated  recurrence  given  by  the  last  of 
the  defining  equations  for  {a;},  which  for  the  ILU  factorization  can  be  rewritten  as 

ay  =  a  —  bd/otj-\  -  ce/a;_„.  (3.7) 

We  attempt  to  gain  some  insight  into  this  recurrence  by  examining  a  simplified  version  of  it.  By 
analogy  with  the  result  for  the  first  block,  let 

a  +  a'2  -  4 (bd  +  ce)  0, 
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which  is  the  larger  root  obtained  from  formally  substituting  the  constant  a  in  place  of  aj.  a;_i 
and  aj-n  in  (3.7)  and  solving  for  a.  Consider  the  alternative  recurrence 

cn  =  a  —  ce/a,  a}  =  a  -  ce/a  -  bd/atj-i .  j  >  2,  (3.9), 

which  is  a  simplification  of  (3.7)  in  which  q;_„  is  replaced  by  a  of  (3.8). 

Theorem  3.2.  Let  a  be  given  by  (3.S).  Then  for  the  values  of  a.  6,  c,  d  and  e  of  either  difference 


scheme,  the  sequence  defined  by  (3.0)  is  convergent  with  limit  a.  The  limiting  values  are 
a  =  2  +  sj2  +  pf  +  p2  for  centeted  differences. 

q  =  2  +  P\+P2  +  C'l  +  (1  +  Pi  +  P2)2  for  u Pwind  differences. 


Proof.  The  values  of  a  for  the  two  schemes  are  obtained  by  substituting  into  (3.8)  the  values  of  a. 
b.  c.  d  and  e  from  Section  2.  The  simplified  recurrence  (3.9)  has  the  same  form  as  the  recurrence 
of  (3.3),  with  a  replaced  by  a  -  ce/a.  For  both  schemes,  a  straightforward  computation  shows  that 
a  -  ce/a  is  greater  than  zero  and  hence  nonzero.  We  apply  Theorem  3.1  in  the  same  manner  as 
above:  if  £  =  —bd/(n  -  ce/a)'  is  greater  than  or  equal  to  —  7.  then  the  sequence  {a;/(u  -  re/a)} 
converges  to  _ 


1  4-  v/l  -  4 bd/(u  -  ce/a)2 


and  { dj  \  converges  to 


(a  -  re/a)  4-  \y'(n  -  ce/a)2  -  4 bit 


(3.10) 


The  condition  £  >  -  q  is  equivalent  to  hd  <  f(u  -  re/a)2.  To  see  that  the  latter  inequality 
holds  for  both  schemes,  first  note  that  a  >  2  4-  v'-  >  9.  For  the  centered  scheme,  we  wish  to  show 
that 

1-PI  <  4  -  J(l-~  P'-]  +  l-  ■  (3.11) 

If  ;>2  >  1.  then  (3.11)  holds  trivially.  If  p>  <  1.  then  0  <  1  -  />.t  <  1.  and  the  right  hand  side  of 
(3.11)  is  greater  than  4  -  )+2  >  3.  from  which  the  inequality  follows.  For  'he  upwind  scheme,  the 

condition  £  >  —  1  is  equivalent  to 

1  /  1  -s 


1  4-  2pi  <  -(44-  2(;>1  -e  J>2  I 

4 


-P-2 


(3.12) 


But  (1  +2p2)/a  <  — =— *  4 -  p-2-  Inequality  (3.12)  'hen  follows  from  direct  computation. 

2  +  v  2 


It  remains  tr 


that  a  =  a.  i.e.  t hat 


(il  -  CP / a)  -*■  y  (,;  —  ce/a)2  —  Abd 


or  equivalently,  that 


/(a  —  ce/ ct)2  —  4bd  =  2a  —  (a  —  ce/a ). 


(3.13) 


Rewriting  the  right  hand  side  of  (3.13)  as  (2 a2  —  aa  +  ce)/a.  it  can  be  shown  by  direct  computation 
that  this  quantity  is  positive  for  both  schemes.  Hence,  it  suffices  to  show  that  the  square  of  both 
sides  of  (3.13)  are  equal,  which  simplifies  to 

a2  -  act  4-  (bd  +  ce)  =  0. 

The  solution  to  this  equation  is  a  of  (3.8),  which  completes  the  proof. 


Note  that  we  made  use  of  three  'nequalities  in  this  proof: 

1.  a  —  ce/a  >  0,  which  ensures  that  £  is  well-defined  and  figures  in  the  choice  of  the  positive 
square  root  of  (a  —  ce/a)2  in  (3.10): 

2.  2a  —  (a  —  ce/a)  >  0,  which  allows  us  to  square  both  sides  of  (3.13): 

3.  £  >  -i,  so  that  Theorem  3.1  can  be  applied. 

For  the  MILU  factorization,  we  rewrite  the  last  recurrence  of  (3.2)  as 

ay  =  a  —  b(d  4-  e)/ay_i  -  c(d  +  e)/ay_„,  (3.14) 

The  roots  of  the  quadratic  equation  obtained  by  substituting  a  for  ay.  a;_i  and  ay_„  in  (3.14)  are 

a±  \/a2  -  4(6  +  c)(d  4-  e) 


- — - - — - — — - (3.15) 

Let  a  denote  the  root  with  larger  modulus.  As  above,  we  examine  the  simpler  recurrence  derived 
by  replacing  )a;_„}  with  a: 


a\  =  a  —  c(d  +  e)/a,  a;  =  (a  —  c(d  ■+■  e)/a)  -  b(d  +  e)/a;_] .  j  >  2. 
We  have  the  following  convergence  result: 


(3.16). 


Theorem  3.3.  For  either  difference  scheme.  let  a  denote  the  larger  root  given  by  (3.15)  and  let 
£  =  —b(d  +  e)/(a  —  c(d+  e)  /a)2.  For  the  values  of  a.  b.  c.  d  and  e  of  the  upwind  scheme,  the  sequence 
(3.16)  is  convergent  with  limit  a.  For  the  values  of  the  centered  difference  scheme,  the  sequence 
(3.14)  is  convergent  to  a  provided  that  £  >  -q.  n  -  c(d  +  e)/a  >  0  and  2a  —  a  +  c(d  +  e)/a  >  0. 


The  limiting  values  are 


a  =  2  +  |pi  +  p2\  for  centered  differences. 

a  =  2(1  4-  pi  4-  />z)  for  upwind  differences. 


0 


► 

I 


Proof.  The  values  of  a  for  the  two  schemes  are  obtained  by  substituting  into  (3.15)  the  values  of  a.  b. 
c,  d  and  e  from  Section  2.  For  the  upwind  scheme,  as  in  the  proof  of  Theorem  3.2,  a  —  c(d  +  e)/a  >  0 
follows  from  direct  substitution.  The  condition  £  >  —  i  is  equivalent  to 


2(1  +  -Pi)  <  -  4  +  2/>i+2p2  - 


2(l  +  2p2) 


a 


-  2  +  p,  +  p2  - 


1  +2p2 


a 


But  a  >  2,  so  that 

1  +  2»2  3 

2  +  Pi  +  P2 - >  «  +  Pi  ■ 

a  1 

Consequently,  it  suffices  to  show  that  (§  +  pi)2  >  2(l-t-2pi),  which  follows  directly.  For  the  centered  ' 

scheme,  we  are  only  concerned  with  values  of  p\  and  p2  for  which  a  —  c(d  +■  e)/a  >  0  and  £  >  —  1. 

Applying  Theorem  3.1  to  both  difference  schemes,  the  sequence  defined  by  (3.16)  is  convergent 

with  limit _ 

fa  -  c(d  +-  e)/a)  +  \J (a  —  c(d  +  e)/a)'2  -  4 b(d  +  e) 
a  ~  9  ' 

The  proof  that  a  —  a  is  identical  to  the  analogous  proof  of  Theorem  3.2.  The  condition  2a  -  a  + 

c(d  +  e)/a  >  0  again  follows  for  the  upwind  scheme  from  direct  computation. 


The  extra  assumptions  made  for  the  centered  difference  scheme  in  Theorem  3.3  are  not  valid 
for  all  real  p\  and  p2-  The  following  result  gives  sufficient  conditions  for  the  inequalities  to  hold. 
We  defer  a  proof  to  the  Appendix.  An  illustration  of  these  values  is  shown  in  Figure  1.  Note 
that  for  most  values  where  Lemma  3.1  does  not  hold,  either  p\  or  p2  is  large,  and  these  values  are 
of  somewhat  limited  practical  interest  (l9j.  See  the  Appendix  for  some  other  comments  on  these 
values. 

Lemma  3.1.  For  the  values  of  u.  b.  c.  d  and  e  of  centered  difference  scheme,  the  inequalities  £  >  -  1. 
a  -  e[d  +  e)/ct  >  0  and  2a  -  a  +  c(d  +  e) /a  >  0  hold  for  all  p\  and  p2  satisfying 

1.  3  <  p2  <  S  and  p\  >  -  1:  or 

2.  -  1  <  p2  <  3  and  p\  arbitrary:  or 

3  P'2  <  -  1  and  either  p\  <  1(  1  —  po)  or  pi  >  2  —  p2. 

We  define  the  constant  coefficient  ILU  and  MILU  factorizations  as  in  (3.1).  except  we  take 
the  diagonal  entries  (a;  )  to  have  the  limiting  values  jo)  of  Theorems  3.2  and  3.3.  We  denote  the 
off-diagonal  entries  J;.  *;y,  and  tj}  of  the  constant  coefficient  factorizations  by  .j.  e  and  t/. 
respectively. 
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p  i 

p  =  Ph 

Eigenvalues  of 
Symmetric  Part 
Leftmost  Rightmost 

Real  Parts  of  Eigenvalues 
of  Operator 

Leftmost  Rightmost 

Iterations  of 
Orthomin(l) 

Iterations  .-.f 

G  MR  ESI  20) 

20 

.625 

.  159E0 

.1 11E1 

.694E0 

.1 10E1 

19 

11 

30 

.9375 

.696E0 

.102E1 

.949E0 

.  102E1 

6 

6 

40 

1.25 

-.148E1 

.548E1 

.939E0' 

.119E1 

17 

8 

50 

1.5625 

-.493E2 

.561E2 

.366E0' 

.138E1 

DC 

11 

GO  ! 

1.375 

-.393E3 

.401E3 

.300E0’ 

.152E1 

DC 

13 

LOO 

2.1375 

-.393E5 

.393E5 

.621E0* 

.219E1 

DC 

27 

L50 

1.6875 

-.520E6 

.520E6 

.466E0* 

.702E1 

DO 

74 

L75 

5.4688 

-.115E7 

.115E7 

-.989E2 

.445E1 

DC  ! 

>  100 

100 

6.25 

-.223E7 

222  E7 

-.775E4 

.553E1 

DC 

>  100 

125 

7.0313 

-.495E7 

.477E7 

-.177E6 

Overflow 

DC 

DC 

Table  3:  Eigenvalues  and  performance  of  iterative  solvers. 


centered  differences,  P\  —  Pi  =  P.  h  =  1/32,  ILU. 

leftmost  eigenvalues  of  the  symmetric  part  change  in  sign  and  the  performance  of  Orthomin(l) 
degrades.  The  diminished  effectiveness  of  GMRES(20)  coincides  more  closely  with  instability  of 
the  preconditioning  than  in  the  previous  example,  although  the  eigenvalues  of  A  again  appear  to 
be  less  sensitive  than  those  of  the  symmetric  part.  * 


p  ! 

1 

p  =  Ph 

Eigenvalues  of 
Symmetric  Part 
Leftmost  Rightmost 

Real  Parts  of  Eigenvalues 
of  Operator 

Leftmost  Rightmost 

Iterations  of 
Orthomin(  1) 

Iterations  of 
GMRES(20) 

50 

1.5625 

.433E-1 

.156E1 

.863E0 

.135E1* 

32 

19 

60 

1.875 

.415E-1 

.177E1 

.307E0* 

.  129E1* 

32 

19 

100 

2.1875 

.340E-1 

.380E1 

.562E0* 

.219E1 

43 

31 

110 

3.4375 

.320E-1 

.512E1 

.575E0' 

.264  El 

14 

13 

L20 

3.75 

.2S9E-1 

.744E1 

.527EO- 

.327E1 

63 

55 

130 

4.0625 

-.17GE1 

.121E2 

.547E0* 

418E1 

S3 

76 

140 

4.375 

-999E1 

.242E2 

.554E0' 

548E1 

DC 

98 

150 

4.6875 

-.481E2 

.681E2 

399E0* 

.702E1 

DC 

>  100 

175 

5.4688 

-.261E4 

.251E4 

-.989E2 

.445E1 

oc 

>  100 

100 

6.25 

-.764E5 

.687E5 

-.775E4 

553E1 

1 

DC 

>  100 

125 

7.0313 

-.125E7 

.107E7 

177EG 

Overflow- 

DC- 

* - i 

>  100 

Table  4:  Eigenvalues  and  performance  of  iterative  solvers. 


centered  differences.  —  P\  =  Pi  =  P.  h  —  1/32,  ILU. 


Tables  5  and  6  show  the  results  for  problems  from  the  third  and  fourth  classes,  resoeetivelv : 
centered  differences,  MILU  preconditioning,  and  either  p\  =  pi  =  p  >  0  or  p  —  —p\  —  pi  >  0.  For 

"Of  course,  the  real  parts  are  not  the  only  indicators  of  large  eigenvalues.  Although  the  Chebvshev- Arnoldi  method  is  not 
specifically  designed  to  fine!  eigenvalues  with  large  imaginary  parts,  the  Arnoldi  computation  docs  compute  a  set  of  estimates 
whose  real  parts  lie  between  the  extreme  real  parts.  We  did  not  observe  anv  such  eigenvalue  estimates  with  imaginary  parts  ) 

that  were  significantly  larger  than  their  extreme  r**al  parts  for  either  of  the  first  two  examples 
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and  IJ~T  v  have  the  same  stability  properties  as  L~lv  and  U~l  v.  respectively.)  The  stopping  cri¬ 
terion  for  the  eigenpair  estimates  (A,  v)  is  ||.4i’  —  Au||  <  10~6  where  ||t’||  =  1.  (For  entries  in  the 
tables  below  marked  with  an  asterisk  (*),  convergence  of  the  eigenvalue  computations  was  slow 
and  the  stopping  criterion  was  not  satisfied.  In  these  cases,  the  extreme  eigenvalues  are  not  well 
separated,  but  the  values  shown  give  a  reasonable  idea  of  the  approximate  values  of  the  set  of 
extreme  eigenvalues  [16].) 

For  positive  integer  k,  Orthomin(k)  and  GMRES(k)  generate  a  sequence  of  approximate  solu¬ 
tions  {-Tj}  to  (2.2)  that  minimize  ||g  —  Ax} |]  over  a  space  of  dimension  at  most  k+1  [6,  17].  Recall 
that  Orthomin(k)  is  known  to  converge  only  when  the  coefficient  matrix  has  positive-definite  sym¬ 
metric  part  [6],  In  our  experience,  Orthomin(k)  is  more  robust  when  more  directions  are  used: 
we  use  k=l  to  try  to  identify  when  a  preconditioning  is  weak.  GMRES(k)  will  solve  arbitrary 
nonsingular  problems  for  large  enough  k,  although  for  any  given  value  of  k  it  is  only  guaranteed  to 
compute  a  sequence  of  iterates  with  nonincreasing  residuals.  For  testing  these  methods,  the  right 
hand  side  of  (2.1)  was  chosen  so  that  (2.4)  is  the  continuous  solution,  and  the  initial  guess  for  the 
discrete  solution  was  zero,  as  in  Section  2.  In  the  examples  below,  we  allowed  these  methods  to 
run  for  at  most  100  iterations,  and  oo  indicates  that  the  residual  norms  {||r,||}  stopped  decreasing 
at  some  point  during  the  run. 

We  first  consider  examples  from  the  four  classes  of  problems  derived  from  centered  differences 
at  the  end  of  Section  4.  In  Table  3,  we  treat  the  first  class:  centered  differences.  p\.  p2  >  0.  ILU 
preconditioning.  For  fixed  h  =  1/32,  we  examine  various  values  of  P  =  P\  —  Pi  >  0.  with  p  =  Ph. 
The  stability  boundary  for  this  set  of  problems  is  p  —  1.  The  most  dramatic  correlation  between 
stability  and  performance  is  in  the  eigenvalues  of  the  symmetric  part.  As  p  increases  through  1. 
the  leftmost  computed  eigenvalues  change  from  positive  to  negative,  and  both  extreme  eigenvalues 
grow  rapidly  as  p  increases.  Failure  of  Orthomin(l)  to  converge  to  the  correct  solution  coincides 
almost  exactly  with  the  set  of  values  p  giving  negative  eigenvalues  for  the  symmetric  part.  The 
eigenvalue  computations  for  the  matrix  ,4Q-’  itself  appear  to  be  less  sensitive  to  stability,  although 
for  p  >>  1  the  leftmost  real  parts  also  become  large  and  negative.  Similarly,  the  performance  of 
GMRES(20)  degrades  for  large  values  of  p.  but  it  is  less  sensitive  than  that  of  Grthomin(l). 

In  Table  4,  we  consider  examples  from  the  second  class  of  problems  of  Section  4:  centered 
differences.  p\  <  0,  pi  >  0.  ILU  preconditioning.  Again,  we  fix  h  =  1/32  and  vary  P  —  —  P\  = 
Pi  >  0.  The  stability  bound  for  this  set  of  problems  is  p  =  Ph  =  2  +  ^  3.732.  The  results 

are  similar  to  those  of  the  previous  example.  As  p  increases  through  the  stability  bound,  the 
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We  return  to  the  question  of  multiple  roots  for  MILU  and  centered  differences.  Recall  that 
the  characteristic  polynomials  may  have  multiple  roots  in  the  second  or  fourth  quadrants  for  some 
choices  of  p\  and  p 2  for  which  both  have  modulus  greater  than  one.  However,  as  we  have  just 
shown,  the  requirement  that  the  maximal  root  be  bounded  by  one  is  also  violated  for  these  values 
of  pi  and  p2 ,  so  that  the  presence  of  multiple  roots  does  not  affect  the  stability  analysis. 

5.  Upwind  differences.  By  assumption,  we  are  examining  the  upwind  scheme  only  for  pi  >  0, 
P2  >  0.  For  the  lower  triangle  of  both  the  ILU  and  MILU  factorizations,  iS  =  -(1  +  2pi).  7  = 
-(1  +  2 P2 ) -  Both  these  quantities  are  negative,  so  that  Case  1  of  Theorem  4.1  applies.  For  ILU, 
q  =  1  +  s  +  v'l  +  s2,  where  £  =  1  +  pi  +  P2  >  1,  and  A„(l)  =  1  —  £  +-  \/ 1  4-  £2  >  1.  Hence,  the  lower 
triangular  solve  is  always  stable.  For  MILU,  a  —  2£,  so  that  A„(l)  =  0  and  the  lower  triangular 
solve  is  also  always  stable.  For  the  upper  triangle  of  both  factorizations,  5  =  —  1  /a,r/  =  —1/a,  so 
that  Case  1  of  Theorem  4.2  applies.  But 


p„(l)=l--  = 
a 


1  - 


>  1  - 


1  +  $+  x/1  +  s2  '  2+  y/2 

9 


>  0 


1 - >0. 

9  c  — 

-s 

Hence,  for  the  upwind  scheme,  both  upper  triangular  solves  are  always  stable. 


for  ILU 
for  MILU 


5.  Correlation  of  Numerical  Performance  and  Stability 

In  this  section,  we  show  that  there  is  a  correlation  between  the  numerical  properties  of  the 
true  preconditioned  operators  and  4e  stability  of  the  constant  coefficient  preconditioning  solves. 
Let  .4  denote  the  preconditioned  operator.  .4Q_1.  For  various  values  of  the  parameters  p\  and  p2 
(determined  by  P\.  P>  and  h).  we  examine  four  properties  of  the  preconditioned  matrix  and  linear 
system: 

1.  the  extreme  eigenvalues  of  the  symmetric  part  (.4  +  AT)/2: 

2.  the  extreme  real  parts  of  the  eigenvalues  of  .4; 

3.  the  performance  of  Orthomin(l)  with  preconditioning  by  the  incomplete  LU  factorizations: 

4.  the  performance  of  GMRES(20)  [17]  with  preconditioning  by  the  incomplete  LU  factorizations. 
All  computations  were  performed  on  a  VAX11-780  in  double  precision  (55  bit  mantissa). 

The  eigenvalues  were  computed  by  Arnoldi’s  method  with  Chebvshev  acceleration,  which  can 
compute  eigenvalues  with  algebraically  largest  or  smallest  real  parts  efficiently  if  these  eigenvalues 
are  well  separated  from  the  interior  eigenvalues  [16].  This  method  repeatedly  computes  matrix- 
vector  products  of  the  form  .4Q-1e  and,  for  the  symmetric  part.  (.4Q-1)r,  so  that  the  precondi¬ 
tioning  triangular  solves  figure  prominently  in  the  computations.  (The  transpose  operations  L~r  r 
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4.  Centered  differences,  pi  <  0,  p2  >  0,  MILU.  For  this  class  of  problems,  7  <  0,  and  either 
3  <  0  for  - 1  <  pi  <  0  and  Case  1  of  Theorem  4.1  applies;  or  3  >  0  for  pi  <  -1  and  Case  3  applies. 
In  the  first  instance, 

r  0  if  P2  >  -pi 

An(l)  =  Qr  +  /3  +  7=  < 

l  “  2(pi  +  p2)  if  P2  <  ~Pl- 

Both  these  expressions  are  nonnegative,  so  that  the  lower  solve  is  stable  for  -1  <  p\  <  0.  In  the 
second  instance, 

{  2(1  +  pi)  if  p2  >  -Pi 

An(-1)  =  a  -  J  +  7  =  < 

1  2(1  -  p2)  if  p2  <  -pi- 

The  first  expession  is  negative  for  p\  <  —  1,  and  the  second  expression  is  nonnegative  for  0  <  P2  <  2. 
Thus,  the  lower  solve  is  stable  for  -1  <  px  <  0  or  0  <  p2  <  1.  The  analysis  for  the  upper  triangle 
is  similar  and  gives  the  same  stability  region. 


-«  -2  0  2  1  8 
Figure  8:  Stability  regions  for  the  MILU  factorization. 

The  stability  regions  for  problem  classes  3  and  4  are  shown  in  the  first  and  second  quadrants 
of  Figure  8.  For  both  ILU  and  MILU.  the  stability  regions  in  the  third  quadrant  are  the  reflections 
over  the  diagonal  line  p2  =  —p\  of  the  regions  in  the  first  quadrant,  and  the  stability  regions  in 
the  fourth  quadrant  are  the  reflections  over  the  line  p2  =  Pi  of  the  regions  of  the  second  quadrant. 
This  can  be  seen  by  replacing  p\  and  p2  by  —71  and  — <72,  respectively  in  the  analysis  of  the  four 
examples  above.  Figures  7  and  8  show  the  full  stability  regions. 


The  stability  region  for  this  problem  class  is  the  intersection  of  the  regions  for  the  lower  part 
and  the  upper  part.  This  region  is  shown  in  the  first  quadrant  in  Figure  7.  The  dashed  curve  is 
the  boundary  of  the  stability  region  for  the  upper  triangle,  showing  that  the  lower  solve  has  more 
stringent  stability  restrictions  than  the  upper  solve.  In  the  simple  case  of  pj  =  p2  =  p,  the  stability 
bound  is  p  <  1. 

2.  Centered  differences,  p\  <  0,  p2  >  0,  ILU: The  analysis  for  both  the  lower  triangle  and  the 
upper  triangle  is  essentially  the  same  as  that  given  for  the  upper  triangle  in  the  previous  example. 


For  both  solves,  the  stability  region  consists  of  the  set  of  points  (pi,P2)  in  the  second  quadrant  of  ’ 
R2  satisfying 

_  ^  2P1  +  1 

P  ‘2  ^  - • 

Pi  +  2 

This  region  is  shown  in  the  second  quadrant  of  Figure  7.  In  the  case  of  —  pl  =  p2  =  p,  the  stability 
bound  is  p  <  2  +  \/3. 


Figure  7:  Stability  regions  for  the  ILU  factorization. 


3.  Centered  differences,  p\,  p 2  >  0.  MILU.  For  the  lower  triangle.  Case  1  of  Theorem  4.1 
applies,  and  A „ ( 1 )  =  0.  Hence,  the  lower  triangular  solve  is  always  stable.  For  the  upper  triangle, 
all  four  cases  can  occur,  and  the  solve  is  stable  in  each  case.  The  analysis  for  all  four  cases  is 
straightforward.  We  present  Case  3  as  representative:  <5  >  0  (pi  >  1)  and  7  <  0  (0  <  p2  <  1).  The 
solve  is  stable  for  even  n  if  and  only  if  pn(— 1)  >  0.  But 

*Hn(- 1)  =  2(1  +p2)  >  2  >  0. 
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For  the  upper  triangle,  all  four  cases  of  Theorem  4.2  can  occur,  with  6  <  0  if  and  only  if  pi  <  1 
and  rj  <  0  if  and  only  if  P2  <  1.  We  examine  Case  2  in  detail,  noting  without  proof  that  the  upper 
triangular  solve  is  stable  for  the  other  three  cases.  Case  2  corresponds  to  {(pi,P2)|pi>P2  >  1}-  By 
Theorem  4.2,  the  solve  is  stable  for  odd  n  when  p„(  — 1)  <  0.  After  scaling  by  a ,  this  is  equivalent 
to 

pi  +  p2  -  4  <  \j2  +  p\+p\.  ( 4.8) 

There  are  now  three  main  subcases.  For  Case  2a,  if  Pi  +  P2  —  4  <  0,  then  (4.8)  is  trivially  true. 
Otherwise,  squaring  both  sides  and  simplifying  gives 


P2 (pi  -  4)  <  4pi  -  7 


as  the  condition  for  stability.  This  is  trivially  true  when  pi  =  4.  Cases  2b  and  2c  are  determined 
by  the  two  branches  of  the  hyperbola  p2  =  (4pi  —  7)/(pi  —  4).  Case  2b  consists  of  the  set  cf  points 
(pi  i  P2)  in  the  upper  right  quadrant  of  R2  for  which  p\  <  4  and  P2  >  (4pi  -  7)/(pi  -  4).  Case  2c 
consists  of  the  points  in  the  upper  right  quadrant  for  which  pi  >  4  and  P2  <  (4p!  -  7)/(pi  -  4). 
A  diagram  of  the  stability  region  in  the  P1-P2  plane  for  the  upper  solve,  with  labels  for  the  various 
cases  and  subcases,  is  given  in  Figure  6.  (The  hyperbola  branch  for  Case  2b  is  not  shown.) 


Figure  6:  Labeled  stability  regions,  pi,  P2  >  0,  upper  triangle, 
ILU. 
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Hence,  with  the  alternative  ordering.  Case  4  is  equivalent  to  Case  3  with  the  roles  of  3  and  7 
interchanged.  The  lower  triangular  solve  is  stable  for  even  n  if  and  only  if  A„(-l)  >  0.  (The  proof 
that  A„  also  has  no  multiple  roots  is  the  same  as  that  for  An;  we  omit  the  details.) 

VVe  summarize  these  observations  as  follows: 

Theorem  4.1.  Necessary  and  sufficient  conditions  for  the  lower  triangular  solve  to  be  stable  are: 

1.  for  3  <  0.  7  <  0:  A„(l)  =  a  +  3  +  7  >  0; 

2.  for  3  >  0,  7  >  0  and  n  odd:  An(— 1)  =  —  a  +  f3  - f  7  <  0; 

3.  for  3  >  0.  7  <  0  and  n  even:  A„(-l)  =  a  -  3  +  7  >  0: 

4.  for  3  <  0,  7  >  0  and  n  even:  An(— 1)  =  a-  ^  +  3>0. 

The  identical  analysis  can  be  used  to  determine  the  stability  of  the  upper  triangular  solve, 
based  on  the  largest  characteristic  root  of  pn  in  (4.5).  We  again  distinguish  among  four  cases, 
depending  on  the  signs  of  <5  and  rj. 

Theorem  4.2.  .Yecessary  and  sufficient  conditions  for  the  upper  triangular  solve  to  be  stable  are: 

1.  for  8  <  0,  p  <  0:  p„(l)  =  1  +  6  +  i)  >  0; 

2.  for  8  >  0,  p  >  0  and  n  odd:  p„(~l)  =  -1  ~h  8  +  r/  <  0; 

3.  for  6  >  0,  rj  <  0  and  n  even:  /r„(-l)  =  1  -  8  +  n  >0; 

4.  for  8  <  0,  r]  >  0  and  n  even:  /i,.(-l)  =  l-  /?+6>0.  where  iin(z)  =  zn  +  r]zn~l  +  8. 

Given  values  of  px  and  P2  (determined  by  Px.  P2  and  h)  and  choice  of  difference  scheme. 
Theorems  4.1  and  4.2  can  be  used  to  determine  whether  the  ILU  or  MILU  factorization  results  in  a 
stable  preconditioner.  We  now  characterize  some  particular  classes  of  problems  and  factorizations. 
As  we  will  show  in  Section  5,  the  conclusions  of  the  Theorems  also  appear  to  hold  for  the  parities  of 
n  not  covered  by  the  analysis.  Hence,  in  the  following  we  do  not  limit  our  conclusions  to  particular 
parities.  Recall  that  for  the  centered  difference  discretization,  3  =  -(1  +  pi),  7  =  -(1  +  P2). 
8  =  (-  I  +  Pi)/a  and  p  =  (- 1  +  p2)/<*- 

1.  Centered  differences ,  p\,  pi  >  0.  ILU.  By  Theorem  3.2,  a  =  2  +  yj2  +  p\  +  p\.  For  the 
lower  triangle,  3  <  0  and  7  <  0,  so  that  Case  1  of  Theorem  4.1  applies.  The  lower  triangular  solve 
is  stable  if  and  only  if 

An(l)  =  \J2  +  pj  +  pi  -  (p\  +  P2)  >  0.  (4.7) 

.After  simplifying,  (4.7)  reduces  to  p\pi  <  1.  which  is  the  necessary  and  sufficient  condition  for 
stability  of  the  lower  triangular  solve. 
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Figure  5:  Shape  of  characteristic  polynomials.  3  >  0.  ,  >  0. 

Left  side:  n  odd.  Right  side:  n  even. 

where  3  and  -7  are  less  than  or  equal  to  0.  Up  to  a  sign,  the  polynomial  on  the  right  of  (4.6)  has 
the  form  considered  for  Case  1.  Hence,  the  analysis  for  Case  1  implies  that  for  odd  n,  the  lower 
triangular  solve  is  stable  if  and  only  if  A„(-l)  <  0. 

Similarly,  for  Case  3,  when  n  is  even, 

A„(-*)  =  otzn  +  {-,8)zn~l  +  -7, 

which  again  has  the  form  studied  for  Case  1.  Thus,  in  Case  3  with  n  even,  the  lower  triangular 
solve  is  stable  if  and  only  if  An(  —  1)  >  0. 

Finally,  for  Case  4,  we  consider  a  recurrence  analogous  to  (4.2)  corresponding  to  a  reordering 
of  unknowns.  Recall  that  in  Section  2  we  assumed  that  the  unknowns  are  ordered  by  lines  along 
the  x-direction.  For  the  lower  triangular  solve  of  (4.1).  the  computations  for  the  unknowns  17  are 
unchanged  (except  for  the  order  in  which  they  are  done)  if  the  unknowns  are  computed  along  the 
lines  in  the  y-direction.  If  we  renumber  the  unknowns  for  the  lower  triangular  solve  according  to 
this  ordering,  then  the  typical  recurrence  in  the  new  ordering  is 

avj  +  -717-1  +  3vj-„  =  wj, 

and  the  corresponding  characteristic  polynomial  is 

Am(-)  =  ac"  +rr"- 1  +3. 
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The  same  argument  works  for  the  upwind  scheme,  for  both  ILU  and  MILU.  For  /in,  the  only 
possiblity  for  a  nonzero  multiple  root  is  z **  =  -(n  -  1  )6/n,  and  the  identical  argument  shows  that  * 

|s**|  <  1.  For  MILU  and  centered  differences,  the  same  reasoning  shows  that  there  are  no  multiple 
roots  if  pi  and  pi  have  the  same  sign,  or  if  either  |pi)  £  1  or  |p2|  <  1.  There  can  be  multiple  roots 
otherwise.  However,  we  will  show  at  the  end  of  this  section  that  such  roots  play  no  role,  so  in  the 
following  we  will  ignore  the  possibility  of  multiple  roots. 

It  is  sufficient,  then,  to  examine  the  moduli  of  the  largest  roots  of  the  two  characteristic 
polynomials.  For  A„,  we  distinguish  among  four  cases: 

1.  3  <  0,  7  <  0;  3.  3  >  0,  7  <  0; 

l 

2.  3  >  0,  7  >  0;  4.  3  <  o,  7  >  0. 

For  Case  1  and  nonzero  3  and  7,  it  can  be  shown  using  elementary  calculus  that  the  graph  of  the 

real  values  of  A„  has  one  of  the  two  shapes  given  in  Figure  5,  depending  on  whether  n  is  odd  or 

even.  In  particular,  An(0)  =  7  <  0,  A„(<)  has  a  local  minimum  at  tm  —  >  0,  and  An(<)  1 

is  strictly  decreasing  for  t  between  0  and  tm  and  strictly  increasing  for  t  >  tm.  Therefore.  A„  has 

precisely  one  positive  real  root,  r,  and  A„(<)  >  0  for  all  t  >  r.  The  same  argument  works  for  3  <  0, 

7  =  0.  If  3  =  0  and  7  <  0,  then  X„(t)  is  strictly  increasing  for  t  >  0  and  again  has  precisely  one 

positive  real  root,  r.  (When  3  —  7  =  0,  0  is  an  n-fold  root  of  An  and  the  solve  is  stable.  In  this  ' 

trivial  case,  the  lower  triangular  matrix  is  actually  a  diagonal  matrix  and  we  will  not  discuss  it 

further.) 

We  claim  that  the  largest  positive  root  r  is  a  root  with  largest  modulus.  For  if  re'9  is  any  root 

l 

with  largest  modulus,  then  by  definition  r  >  r  and 

arne'ne  =  -3rn-lei{n~1)e  -  7. 

Taking  the  modulus  of  both  sides  and  applying  the  triangle  inequality, 

afn  =  |  -  -  7|  <  l/ilf""1  +  |7|, 

i.e.  A„(r)  <  0.  Since  r  >  0,  it  follows  that  r  <  r.  Consequently,  r  =  r.  Whether  or  not  r  is 
greater  than  I  can  be  determined  from  the  sign  of  A„(l).  For  if  A„(l)  <  0,  then  since  1  is  a  positive 
number,  it  follows  that  1  must  lie  to  the  left  of  r,  i.e.  r  >  I.  Conversely,  if  A„(l)  >  0,  then  r  <  I. 

Hence,  for  Case  1,  the  lower  triangular  solve  is  stable  if  and  only  if  An(l)  >  0. 

For  Cases  2  —  4,  we  can  only  give  a  partial  analysis,  which  characterizes  stability  only  for 
certain  parities  of  n.  Consider  Case  2  and  n  odd.  Then  ^ 

An(  — -)  =  -  3zn~l  -  7)  =  -(«"  +  3zn +  7),  (4.6) 
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Consider  the  homogeneous  case,  w  =  0.  If  A„  has  n  distinct  roots  {z\, . . . ,  zn},  then  it  is  well  known 
(see  e.g.  [13])  that  for  j  >  n  the  homogeneous  solution  of  (4.2)  is 

Vj=cxz{-\ - +cnzJn,  (4.3) 

where  cj, . . .  ,c„  are  determined  by  the  initial  values  i/i,. . . ,  vn-  If  some  root  z,  of  An  has  modulus 
greater  than  1,  then  zi  will  grow  as  j  increases,  and  any  error  in  cs  (caused,  say,  by  an  error  in  the 
initial  conditions)  will  result  in  increasing  errors  in  the  solution  (4.3).  If  a  root  z,  has  multiplicity 
m,  then  the  homogeneous  solution  also  contains  linear  combinations  of 

J'0-I)-  (j -t+l)zi~l,  t  —  1, . . . ,  m  —  1  (4.4)' 

as  components;  if  z,  has  modulus  greater  than  or  equal  to  1,  then  any  errors  in  the  coefficients  of 
(4.4)  will  also  be  enhanced  with  increasing  j.  Therefore,  we  say  that  the  lower  triangular  solve  is 
stable  if  all  the  roots  of  its  characteristic  polynomial  are  less  than  or  equal  to  1  in  modulus  and  no 
root  with  unit  modulus  has  multiplicity  greater  than  1.  It  is  unstable  otherwise. 

The  typical  computation  in  the  upper  triangular  solve  is 

Vj  +  bv1  +  {  +  'yvj+n  =  Wj , 

where  v]  +  l  and  t>;+„  are  given.  To  fit  this  computation  into  the  setting  of  recurrences,  we  renumber 
th«  unknowns  so  that  they  are  computed  in  order  of  increasing  indices,  i.e.  let  bj  =  v^-j.  Then 
the  upper  triangular  computation  has  the  form 

Vj  +  bVj-X  +  7^-n  =  w,V-y, 

and  the  associated  <  haiacteristic  polynomial  is 

»„{;)  =  zn+6:n~l  +  n  =  0.  (4.5) 

We  say  that  the  upper  triangular  solve  is  stable  if  all  the  roots  of  its  characteristic  polynomial  have 
modulus  less  than  or  equal  to  1  and  no  root  with  unit  modulus  has  multiplicity  greater  than  1.  It 
is  unstable  otherwise. 

We  first  note  that  for  the  ILU  factorization  and  for  the  MILU  factorization  with  the  upwind 
scheme,  there  are  no  multiple  roots  of  unit  modulus.  Any  multiple  root  of  A n(z)  must  also  be  a 
root  of  A'n(a)  =  zn~2{naz  +  (n  -  l)/tf).  Thus,  the  only  possibility  for  a  nonzero  multiple  root  is 
z *  =  — (n  -  l)[3/na.  For  ILU  and  the  centered  difference  scheme, 

i**i<!£!  =  — 1.1  +  Plj  ■,■<!■ 

«  2+  s/2  +  p\  +  pj 
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We  have  not  seen  any  examples  in  which  the  sequences  (3.9)  and  (3.16)  are  convergent  but  the 
true  diagonal  sequences  {ay}  are  not  convergent.  However,  we  do  have  cases  where  convergence 
of  the  diagonal  sequence  is  slower  than  in  the  examples  above.  For  example,  Figure  4  graphs  the 
computed  and  limiting  values  for  P[  =  —50,  P2  =  50,  h  =  1/32,  centered  differences  with  the 
MILU  factorization,  from  four  blocks  of  the  factorization.  Within  each  block,  the  diagonal  values 
approach  a  limit,  and  as  the  factorization  proceeds,  the  limit  for  each  block  gets  closer  to  the  limit 
of  Theorem  3.3.  However,  the  initial  values  within  each  block  oscillate,  and  these  oscillations  are 
both  larger  in  magnitude  and  occur  at  higher  indices  (mod  n)  within  the  block  at  the  later  stages 
of  the  factorization. 

4.  Stability  Analysis 

In  this  section,  we  examine  the  numerical  stability  of  the  lower  triangular  and  upper  triangular 
solves  performed  with  the  constant  coefficient  factorizations  introduced  in  Section  3.  Given  an 
incomplete  factorization  LU ,  the  preconditioning  operation  consists  of  a  pair  of  triangular  solves 
of  the  form 

Lv  =  w,  Uv  =  w.  (4.1) 

Typically  these  operations  are  performed  once  each  per  iteration  of  the  basic  iterative  method. 
Consider  the  lower  triangular  solve.  The  typical  computation  for  the  y  'th  entry  of  v  has  the 

form 

vj  =  ~(w}  ~  dvj. ,  -  7V,-")' 

where  i’,_i  and  ty _ „  have  been  previously  computed.  Equivalently,  most  entries  of  v  satisfy  the 
n’th  order  inhomogeneous  recurrence  relation 

otVj  +  3vj_  1  +  7  Vj-„  =  Wj.  (4.2) 

Our  stability  analysis  is  based  on  an  analysis  of  this  recurrence,  whose  solution  is  uniquely  defined 

once  n  initial  values,  say  17 . e„,  are  given.  (We  note  that  as  in  computation  of  the  factors,  not 

every  step  of  the  backsolve  satisfies  this  recurrence,  since  there  are  2n  —  1  cases  corresponding  to  grid 
points  next  to  the  boundary  where  the  computation  of  1  y  is  simpler.)  We  wiil  define  the  stability 
of  the  lower  triangular  solve  in  terms  of  properties  of  the  characteristic  polynomial  associated  with 
(4.2),  which  is  given  by 
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Figure  4:  Convergence  of  diagonal  values.  P\  =  -50,  P>  = 
50.  /i  =  1/32,  centered  differences,  MILU  factorization. 
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Figure  1:  Values  of  p \  and  P2  where  Theorem  3.3  holds  for 
centered  differences  and  MILU. 

The  results  of  the  theorems  do  not  rigorously  prove  that  the  true  factors  converge  in  any  sense 
to  the  constant  coefficient  factors.  The  simplifying  assumptions  that  the  diagonal  values  from  the 
previous  blocks  are  constant  and  that  they  take  on  the  values  derived  from  solving  the  specified 
quadratic  equations  are  not  true.  (For  example,  the  limiting  values  for  the  first  block  are  in  general 
different  from  the  quadratic  roots  used  in  the  theorems.)  Moreover,  the  convergence  results  do 
not  say  anything  about  how  close  the  first  n  values  are  to  the  limits,  which  is  the  only  useful 
information  in  this  context.  Nevertheless,  our  numerical  experience  supports  the  idea  that  the 
constant  coefficient  factors  are  reasonable  approximations.  In  Figure  2.  we  graph  the  computed 
values  of  the  ILU  diagonals  (a; }  and  the  limiting  value  a  of  Theorem  3.2,  for  a  15  x  15  grid 
( h  =  1/16)  and  Pi  =  Pi  =  25  (so  that  pi  =  P2  —  1.5625).  As  expected,  the  limiting  value  for  the 
first  block  is  different  from  a.  but  starting  from  the  third  block,  most  of  the  subsequent  diagonal 
values  are  virtually  indistinguishable  from  a.  Figure  3  graphs  the  analogous  data  for  the  same 
problem  and  difference  scheme  with  the  MILU  factorization.  In  this  case,  the  values  from  the  last 
block  and  the  last  values  within  the  other  blocks  differ  from  the  limit  because  they  satisfy  a  different 
recurrence  (see  (3.2)).  For  the  same  problem  with  the  upwind  schemes,  both  factorizations  show 
the  same  qualitative  behavior.  Similarly,  for  centered  differences  on  a  finer  grid  with  the  same 
values  of  p\  and  p2  (h  —  1/32  and  P[  =  P2  =  50),  the  qualitative  behavior  is  identical. 


all  the  values  considered,  Theorem  3.3  holds.  The  analysis  of  Section  4  predicts  that  there  are  no 
stability  restrictions  for  the  problems  of  Table  5.  The  numerical  results  are  consistent  with  this: 
all  of  the  extreme  eigenvalues  vary  smoothly  with  p,  and  neither  iterative  method  has  difficulty 
solving  the  preconditioned  problem.  (We  will  comment  below  on  the  change  in  sign  of  the  leftmost 
eigenvalue  of  the  symmetric  part  at  p  =  6.25.)  The  stability  bound  for  the  problems  of  Table  6 
is  p  =  1.  For  this  set  of  problems,  all  four  performance  indicators  change  dramatically  when  p 
increases  through  1.  (Note  in  particular  that  the  values  of  p  are  significantly  higher  in  Table  5  than 
in  Table  6.) 


p 

p  =  Ph 

Eigenvalues  of 
Symmetric  Part 
Leftmost  Rightmost 

Real  Parts  of  Eigenvalues 
of  Operator 

Leftmost  Rightmost 

Iterations  of 
Orthomin(l) 

Iterations  of 
GMRES(20) 

30 

.9375 

.996E0 

BE 

.100E1 

.104E1 

4 

4  1 

50 

1.5625 

.681E0 

Hi 

.780E0 

.100E1 

7 

7 

o 

o 

•H 

2.1875 

.370E0 

.121E1 

.479E0* 

.100E1 

12 

12  ! 

150 

4.6875 

.141E0 

.136E1 

.352E0 

16 

15  | 

200 

6.25 

-.100E-1 

.146E1 

.251E0 

19 

18 

225 

7.0313 

-.663E-1 

.150E1 

.215E0 

20 

19  I 

Table  5:  Eigenvalues  and  performance  of  iterative  solvers. 


centered  differences,  P\  =  Pi  —  P.  h  =  1/32,  MILU. 


a 

m 

Eigenvalues  of 
Symmetric  Part 
Leftmost  Rightmost 

Real  Parts  of  Eigenvalues 
of  Operator 

Leftmost  Rightmost 

Iterations  of 
Orthomin(l) 

Iterations  of 
GMRES(20) 

30 

.9375 

■ 

.323E1 

52 

35 

31 

.9688 

52 

35 

32 

1 

.819E0* 

51 

36 

33 

1.0313 

-.134E4 

.159E3 

.874E1 

oo 

55 

34 

1.0625 

-.138E4 

.138E4 

DO 

>  100 

36 

1.125 

-.287E3 

.299E3 

oo 

>  100 

1.5625 

-.290E11 

-.343E2 

oo 

oc 

Table  6:  Eigenvalues  and  performance  of  iterative  solvers. 


centered  differences.  —P\  =  Pi  =  P.  h  =  1/32,  MILU. 


As  we  noted  in  Section  4,  some  of  the  stability  analysis  of  Theorems  4.1  and  4.2  applies  only  for 
certain  parities  of  n,  but  the  numerical  performance  seems  independent  of  parity.  As  an  example, 
consider  the  case  of  ILU  preconditioning  with  centered  differences  in  the  second  quadrant  (see 
Figure  7).  The  analysis  at  the  stability  boundary  makes  use  of  Case  3  of  Theorem  4.1  and  Case  4 
of  Theorem  4.2,  both  of  which  require  n  to  be  even.  The  performance  of  ILU  shown  in  Table  4  is 
for  n  =  31,  i.e.  n  odd,  suggesting  that  parity  plays  no  role.  Further  evidence  is  given  in  Table  7, 
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p 

p  =  Ph 

Eigenvalues  of 
Symmetric  Part 
Leftmost  Rightmost 

Real  Parts  of  Eigenvalues 
of  Operator 

Leftmost  Rightmost 

Iterations  of 
Orthomin(  1) 

Iterations  of 
GMRES(20) 

maim i 

.309E-1 

.468E1 

.582EO* 

.249E1 

48 

38 

.287E-1 

.659E1 

.555E0* 

58 

55 

130 

BBSS 

.590E0* 

•  Y  ;  'Y;-  . 

86 

75 

140 

-.591E1 

BUB 

.631EO* 

1 

>  100 

81  i 

Table  7:  Eigenvalues  and  performance  of  iterative  solvers,  • 

centered  differences.  -Pi  —  P-2  =  P.  h  =  1/33,  ILU. 


which  shows  that  for  values  of  p  near  the  stability  boundary  the  performance  for  n  =  32  (h=l/33)  . 
is  essentially  the  same.  0 

Our  previous  experience  with  these  preconditionings  has  been  on  problems  with  just  one  first 
derivative  term  present,  for  which  the  coefficient  is  positive.  Both  preconditioners  have  been  very 
successful  in  solving  such  problems  (7,  8],  The  stability  analysis  of  Section  4  suggests  that  neither 
preconditioning  suffers  from  instability  in  these  cases,  or  for  problems  where  upwind  differences  are  • 

used.  Tables  8  and  9  shows  that  the  extreme  eigenvalues  of  the  symmetric  parts  for  some  problems 
of  these  types  are  well  behaved,  as  expected. 


P2 

Pi  =  Pih 

ILU 

Leftmost  Rightmost 

MILU 

Leftmost  Rightmost 

30 

.9375 

.858E-1 

.123E0 

.242E0 

.305E1 

50 

1.5625 

.172E0 

.132E1 

.240E0 

.259E1 

100 

2.1875 

.386E1 

.147E1* 

.322E0 

213E1* 

150 

4.6875 

.437E0 

.152E1* 

398E0 

.192E1* 

200 

6.25 

.477E0 

.153E1* 

455Et-  • 

.178E1 

225 

7.0313 

.494E0* 

.152E1 

.477E0 

.173E1* 

Table  8:  Eigenvalues  of  the  symmetric  part,  centered  differ¬ 
ences.  Pi  =  0,  h  =  1/32. 


P 

p  =  Ph 

ILU 

Leftmost  Rightmost 

MILU 

Leftmost  Rightmost 

30 

.9375 

.553E-1 

.118E1 

.629E0 

.254E1 

50 

1.5625 

.803E-1 

.115E1 

.846E0 

.198E1 

100 

2.1875 

.150E0 

.lllEl 

.951E0 

.  151E1 

150 

4.6875 

.21SE0 

.109E1 

.973E0 

.135E1 

200 

6.25 

.281E0 

.107E1 

.980E0 

126E1 

225 

7.0313 

.309E0 

.106E1 

.983E0 

.124E1 

Table  9:  Eigenvalues  of  the  symmetric  part,  upwind  differ¬ 


ences,  Pi  —  Pj  =  P  >  0.  h  —  1/32. 
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Note  that  we  are  not  addressing  the  question  of  accuracy  of  the  incomplete  factorization,  and 
we  cannot  be  certain  that  it  is  instability  rather  than  inaccuracy  that  is  causing  the  difficulties 
exhibited  in  Tables  3,  4  and  6.  However,  only  in  cases  in  which  at  least  one  of  the  triangular  solves 
is  unstable  in  the  sense  of  Section  4  do  the  computed  eigenvalues  increase  dramatically  as  they 
do  in  these  three  tables.  We  have  encountered  problems  with  unfavorable  eigenvalue  distributions 
that  appear  to  be  due  to  inaccuracy  of  the  incomplete  factorization.  One  such  case  is  the  example 
of  Table  5:  the  eigenvalues  of  the  symmetric  part  turn  negative  when  p  >  6.25,  but  they  do  not 
change  dramatically  in  magnitude  with  small  changes  in  p.  Another  example  is  as  follows:  we 
let  h  =  1/64  and  P 2  =  100  be  fixed  so  that  P2  =  1.5625,  we  vary  P\  <  0,  and  we  use  centered 
differences  and  the  MILU  factorization.  (This  problem  is  a  member  of  the  fourth  class  of  problems 
analyzed  at  the  end  of  Section  4.)  The  extreme  eigenvalues  of  the  symmetric  part  are  shown  in 
Table  10.  The  leftmost  eigenvalues  are  negative  for  all  values  of  p\  considered,  but  they  are  fairly 
well-behaved  until  the  stability  bound  of  p\  =  —  1  is  reached,  after  which  they  quickly  diverge. 


Pi 

Pi  =  P\h 

Leftmost 

Rightmost 

-10 

-.1563 

-.633E0 

-20 

-.3125 

-.992E0 

.509E1 

-40 

-.625 

-60 

-.9375 

-64 

-1.0 

-68 

-1.0625 

-72 

-1.125 

-.715E6 

.716E6 

-76 

-1.1875 

-.152E9 

.152E9 

Table  10:  Eigenvalues  of  the  symmetric  part,  centered  differ¬ 
ences,  P'2  =  100,  h  =  1/64,  MILU. 


We  conclude  with  two  observations  that  we  have  been  unable  to  explain.  First,  in  Tables  3  and 
4,  the  extreme  real  parts  of  the  eigenvalues  of  the  operator  AQ-1  are  the  same  for  large  p,  although 
the  signs  of  p  1  are  opposite.  Indeed,  we  have  observed  in  some  other  tests  with  centered  differences 
and  ILU  that  for  large  pj  and  P2,  the  extreme  eigenvalues  of  AQ~l  appear  to  be  independent  of 
the  signs  of  pi  and  p2 .*  Second,  for  very  large  values  of  the  parameters  pi  and  p2  in  Tables  3.  6 
ind  10,  the  computed  extreme  eigenvalues  of  the  symmetric  part  are  opposite  in  sign  but  nearly 
e^ual  in  magnitude.  We  do  not  have  good  explanations  for  these  phenomena. 

*We  can  show  that  if  |pj|  =  |pj|,  \pi\  =  \pfn\,  then  the  error  matrices  R  =  Q  —  A  and  fl'  s  Q'  -  A*  for  the  constant 


coefficient  factorizations  are  similar  under  a  diagonal  similarity  transformation.  This  result  holds  for  both  the  ILU  and  MILU 
factorizations,  but  we  have  not  been  able  to  make  use  of  this  observation. 


6.  Appendix 


In  this  section  we  prove  Lemma  3.1.  To  simplify  notation,  in  the  proof  we  use  the  symbols  “x” 
and  *y”  instead  of  “pi”  and  “p2.”  The  three  inequalities  to  be  considered  are  then: 

5>-i:  +  + 


a-c(J  +  *)/oO:  4  -  >  Q; 

v  "  2  +-|ar  -r  y| 


2a  -  a  +  c(d  +  e)/a  >  0  :  2| x  +  y)  +  —  +  |" — -J ±  1)1  >  q, 

2+|x+y| 

We  partition  the  plane  into  five  regions: 


(6.1) 

(6.2) 

(6.3) 


1  ■  *  +  y  <  0; 

2.  0<x  +  y<2  and  1  +  y  >  0; 

3.  0<x+y<2  and  1  +  y  <  0; 

4.  2  <  x  +  y  and  1  +  y  >  0; 

5.  2  <  x  +  y  and  1  +  y  <  0. 

These  regions  are  depicted  in  Figures  9  -  11.  Note  that  in  Region  1,  (2  -  (x  +  y))/(2+  |x+  y|)  =  1. 
so  that  all  three  inequalities  are  simpler  in  this  case. 

In  Region  1,  (6.1)  is  equivalent  to 


<t>(x,  y) 


1 

4 


4-  xy  + 


>  0. 


It  is  straightforward  to  show  that  <t>  =  0  and  V<i>  =  0  on  the  line  y  =  1  —  2x,  and  V2c>  is  positive 
semi-definite  everywhere.  But  <j>  is  quadratic,  so  for  any  ,V  =  (x,  y)T  and  A'o  =  (xo,i/o)r  such  that 
yo  =  1  -  2x0, 

0(.Y)  =  l-(X  -  .Y0)T(V2^)(.Y  -  .Yo)  >  0. 

In  Region  2,  the  right  hand  side  of  (6.1)  is  greater  than 


4  -  (1  +  y)(2  -  (x  +  y)). 


which  can  be  seen  to  be  greater  than  or  equal  to  (1  +  x)(2  —  (x  +  y))  by  direct  computation.  In 
Region  4.  the  left  hand  side  of  (6.1)  is  negative  when  x  >  -1,  so  that  the  inequality  holds  trivially. 
The  same  argument  works  for  (all  of)  Region  5.  We  have  not  been  able  to  establish  (6.1)  for  Region 
3  or  for  Region  4  with  x  <  -1,  although  numerical  tests  suggest  that  it  also  holds  in  these  regions. 


Figure  9:  Values  of  x  and  y  where  inequality  (6.1)  holds. 

The  shaded  area  in  Figure  9  shows  the  set  of  points  in  the  plane  for  which  we  have  demonstrated 
that  (6.1)  holds. 

Inequality  (6.2)  holds  immediately  in  Regions  3  and  4,  since  the  quantity  subtracted  is  negative. 
The  inequality  reduces  to  y  <  3  in  Region  1,  and  for  Region  2,  using  the  fact  that  (2—  (x+y))/(2  + 
\x  +  y|)  <  1,  the  same  condition  is  sufficient.  For  Region  5,  after  being  scaled  by  2  +  x  +  y,  the 
right  hand  side  of  (6.2)  can  be  shown  to  be  greater  than 


8  +  3x  +  ;/=8  +  x+  2/  +  2x>  10  4-  2x  >  13, 


since  x  -1-  y  >  2.  —y  >  1  and  x  >  2  —  y  >  3.  Figure  10  shows  the  set  of  points  where  (6.2)  holds. 

Inequality  (6.3)  holds  trivially  in  Regions  2  and  5.  In  Region  1  the  inequality  reduces  to 
y  <  1  -  2x.  In  Region  3,  the  fact  that  (2  -  {x  +  y))/(2  +  |x  +  y|)  <  1  makes  the  right  hand  side 
greater  than  2x  4-  3 y  +  1,  which  is  greater  than  0  for  y  >  — 1(1  +  2x).  This  determines  a  small 
subset  of  the  region.  Finally,  in  Region  4,  2  4-  \x  +  y|  >  4,  so  that  if  —  1  <  y  <  7,  then 


(1  +  y) 


<  2. 


2  +  |x  +  2/| 

Consequently  the  right  hand  side  of  (6.3)  is  greater  than 


2(x  +  y)  +  2(2  -  (x  +  t/))  =  4  >  0. 


Hence,  (6.3)  holds  in  the  shaded  region  of  Figure  11. 


be  convergent.  For  example,  the  values  pi  =  —4  and  pi  =  4  violate  condition  (6.2),  and  with  these 
values  (from  h  =  1/32,  Pi  =  -128  and  P-i  =  128),  the  MILU  diagonal  sequence  appears  not  to 
be  convergent.  Similarly,  the  values  p\  =  2  and  pi  =  —2  violate  condition  (6.3),  and  the  MILU 
sequence  is  also  not  convergent  for  them  (from  h  =  1/32,  Pi  =  64  and  Pi  =  —64). 
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